HBV Model

In this exercise, we will pull together many aspects of what we have learned so far to calibrate the HBV model to weighted data from the Fraser Experimental Forest in Colorado using precipitation/runoff and potential evapotranspiration (PET) data derived from the methods explored in the snowmelt and evapotranspiration modules.

HBV is: In the previous labs, we have worked with single-process models. HBV is more complex process-based model using many parameters to describe the physical processes of the hydrologic cycle https://rpubs.com/sdhakal/848236

This schematic was reproduced from: Durga Lal Shrestha & Dimitri P. Solomatine (2008) Data‐driven approaches for estimating uncertainty in rainfall‐runoff modelling, International Journal of River Basin Management, 6:2, 109-122, DOI: 10.1080/15715124.2008.9635341

HBV versions with a GUI exists, but the data manipulation is more convenient in a coding environment, and we want you to be able to see ‘under the hood’. Instead of dealing with the inner workings of the model too much, we will put the focus on importing and organizing the data, getting the model running, manual optimization.

The HBV realization we will be using is contained in the R package TUWmodel. Please search that package in your browser and familiarize yourself with it (that will be the first step). TUW simply means “Technische Universitaet Wien”, or Vienna University of Technology. The model is slightly different from the original HBV but still similar enough to call it HBV. We also need to understand the arguments and parmeters required for the model.

Type ‘?TUWmodel’ in your console. Note that the package comes with an example dataset, and you can run examples either by copying and pasting the data(), TUWmodel() and plots into your console, or you can find and click ‘Run examples’ under the Examples header and view the script and plot outputs in the Help window. The first example shows how the ‘TUWmodel’ can be run given specifications for a long set of parameters (param). You can provide the model with precip, air temperature, area and ET along with the parameter set and it will simulate SWE, snowmelt and discharge. We can obtain measurements for each of these variables, so we will import all of this data and test the fit of the simulated data.

https://cran.r-project.org/web/packages/TUWmodel/TUWmodel.pdf

Import data:

We will generate a conceptual runoff model starting with:
- RainMeltWeighted_mm –> liquid input to the ground, combination of rain and melt
- PWeighted_mm –> Precipitation measured at the SNOTEL sites, both liquid and frozen. This is NOT water that went into the ground at that point. This is, however, the P input time series for the HBV model.
- SweWeighted_mm –> Snow water equivalent at the SNOTEL sites. We can use this to evaluate our snow routine.
- TempWeighted –> Either daily mean, min, or max. We will use the mean. - Discharge_mm –> Q observed at the outlet.
- ET –> Potential Evapotranspiration calculated using relative humidity and temperature.

Repo link

Download the repo for this lab HERE

Labwork (20 pnts)

Our first steps will be to import and format the data. As you might recall from the previous modules, this step can require the most script and ‘work’, but is critical to a valid model output. In this assignment, we have simplified much of the data collection for you, as the steps are ones you have already completed in previous modules. In the .csv imported below, we started by collecting the same SNOTEL data that you used in the snowmelt module, but downloaded for a greater date range (water years 2018-2022) (date, SWE depth in mm, and daily precipitation (mm)) However, you will notice that the column names of the provided table contain ‘weighted_’.

#rm(list = ls()) # this clears out the Environment ("fresh start!")
#indata <- read_csv("FEF_lumped_daily_2020_2023.csv")

indata <- read.csv("HBV_data/weighted_data_fool_HBV.csv")%>%
  mutate(date = ymd(date)) %>%
  rename(Q_m3_s = m3_s, Q_mm_day = mm_day, runoff_input.mm = input.mm )

str(indata)

'data.frame':   1826 obs. of  17 variables:
 $ date              : Date, format: "2017-10-01" ...
 $ wtr_yr            : int  2018 2018 2018 2018 2018 2018 2018 2018 2018 2018 ...
 $ weighted_swe.mm   : num  13.3 46.7 53.4 53.4 51.2 ...
 $ weighted_precip.mm: num  26.78 7.65 1.91 0 0 ...
 $ weighted_pcumul.mm: num  26.8 34.4 36.3 36.3 36.3 ...
 $ Q_m3_s            : num  0.00862 0.00743 0.00715 0.00704 0.00704 ...
 $ Q_mm_day          : num  0.282 0.243 0.234 0.23 0.231 ...
 $ Tmax_c            : num  5.55 1.95 5.25 12.75 12.75 ...
 $ Tmean_c           : num  0.8 -1.6 -0.7 6.2 6.65 0.25 3.2 2.5 -7.8 -1.8 ...
 $ Tmin_c            : num  -3.95 -5.15 -6.65 -0.35 0.55 ...
 $ RHmin             : num  43.1 60.3 41.2 28.9 27.4 ...
 $ RHmax             : num  90.6 91.8 91.2 86 78.6 ...
 $ month             : int  10 10 10 10 10 10 10 10 10 10 ...
 $ day               : int  1 2 3 4 5 6 7 8 9 10 ...
 $ SWEdiff.mm        : num  0 33.38 6.75 0 -2.28 ...
 $ Pdiff.mm          : num  0 7.65 1.91 0 0 ...
 $ runoff_input.mm   : num  0 0 0 0 2.28 ...

As we noted in the Snowmelt module, factors such as elevation and vegetation will affect the of snow-to-runoff rate. Similarly, temperature and precipitation differ with elevation, meaning that SNOTEL data from a single location may not fully represent conditions across the entire watershed.

To account for these variations, the “weighted” values for this exercise have been adjusted using linear scaling relationships that estimate average conditions across the watershed. The adjustments use elevation-precipitation and elevation-temperature relationships derived from inverse distance weighting to better approximate spatially distributed hydrological inputs rather than relying on a single point measurement.

/precip_spatial_interp.png

If you are interested, you can start exploring spatial interpolation with: 1. Thiessen polygons 2. Inverse distance weighting 3. Kringing methods to determine if any of these are applicable to your study area.

Other columns in this data include: Q_m3_s and Q_mm_day - Discharge (Q) collected at the Fool Creek outlet by USFS, in units of cubic meters per second and mm per day from April until October. Tmax_c, Tmin_c and Tmean_c - Typically, we could find daily mean, max and minimum temperatures in the SNOTEL datasets, however, this particular station is missing temperature data (due to restrictions on technical access in 2020), so we retrieved temperature data from GridMET through Climate Engine. This highlights the importance of evaluating each variable for completeness. It will save you the headache of running the entire workflow, only to find that model outputs cannot be simulated for the later half of 2020 due to missing input data. RHmin, RHmax - Relative humidity daily min and max, also from GridMET, accessed through Climate Engine SWEdiff.mm, Pdiff.mm and runoff_input.mm - all daily outputs of a temperature based snowmelt model (the same as in the snowmelt module). runoff_input is the estimated daily input to the stream from melted snow and liquid precipitation combined.

Q1 (4 pnts) Generate script (including plots) that can be used to evaluate the structure of the dataset and the duration and consistency of the data. The submitted plot should show at least two time series of ‘indata’ (2 variables on y, date on x), but be sure to check every column for yourself. Do you see any potential issues like missing data or unexpected values? ANSWER:

Keep in mind that you can plot more than two variables in a single plot, but they will be most helpful if you group variables with a similar y-scale. For example, cumulative precipitation or SWE values will have a different range than variables that represent daily measurements like ‘input_mm’. Alternatively, you can add a secondary y-axis to your plot.

BLANK +
  geom_point(aes(x = date, y = runoff_input.mm, color = "runoff_input.mm"), size = 0.5) +  # Assign a label for legend
  geom_point(aes(x = date, y = Q_mm_day, color = "mm/day"), size = 0.5) +  # Assign a different label
  geom_point(aes(x = date, y = weighted_precip.mm, color = "weighted_precip.mm"), size = 0.5) +  # Assign a different label
  BLANK(aes(x = date, y = Pdiff.mm, color = "Pdiff.mm"), size = 0.5) +
  scale_color_manual(values = c("runoff_input.mm" = "blue", "mm/day" = "red", "weighted_precip.mm" = "green", "Pdiff.mm" = 'purple')) +  # Customize colors
  labs(color = "Discharge Units") +  # Legend title
  theme_minimal()

Error: object 'BLANK' not found

2,930m. watershed_area_m2 <- 2640000

Note that while the stream is snow-covered, there are no stage/flow measurements being made at this site. For many of our calculations to work, we will not want NA in our data frames. In this data set, if we view the tabular data, the end discharge reading (Q_mm_day), is very similar to the first in the following calendar year. For this example, we will then fill the NA values with the mean of the final and first readings for each winter (Oct - April) NA string.

Let’s try our quick check again:

Next, we need to add daily potential evapotranspiration (PET) to the dataset. The evapotranspiration module covered some of the numerous pre-built functions available in various R packages designed for ET estimation. For instance, the ‘Evapotranspiration’ package provides ET estimates derived from approximately 20 distinct equations or variations. These functions require precise data formatting to ensure compatibility with the function arguments. To help you structure your function inputs similarly, many packages come with example datasets in their documentation. SPEI is another package that offers ET functions. The github repository for this package has been updated fairly recently, which can be important to verify.

We chose to write our own function here so you could ‘see under the hood’.
Not all packages are equally maintained, as they are often developed by researchers or modelers to improve the repeatably of their work, and contribute to the broader scientific community. Once funding ceases, or if the original developer moves on to new projects, a custom package may no longer receive updates or support. As R, RStudio and other supporting packages are updated, a package may become depreciated. If the package still functions correctly in your current R version and dependencies, there’s no immediate reason to stop using it. However, if you are concerned that future versions of R or dependencies might break the packages you use, you can check the development and maintenance history of custom or specialized packages (or write a package or function for yourself as exemplified below!).

Here is our Hargreaves function, adapted from the ‘Evapotranspiration’ package to minimize re-formatting of our dataframe.

Now we will format the inputs so the data is easily read by the function, and run the function.

Now we’ll add daily ET into our original dataframe:

Q2. (4 pts) Generate a summary dataframe that summarizes columns of ‘indata’ by water year and helps you evaluate the water balance for each year. (I recommend summarizing SWE, precipitation, runoff input, discharge, and PET) Make sure your summary makes sense. For example, some variables are most useful when summarized as a mean, while others are more helpful when viewed as an annual sum, minimum or maximum. If you are unsure of where to start, try using using group_by() and summarize().

Q3 (3 pts) The provided PET data was calculated just using min and max daily air temperature using a Hargreaves approach. View your annual totals for PET and compare it to the amount of precipitation. Are those values reasonable? What can you say about the timing of water demand (PET) and availability (melt and rain)? IMPORTANT: Use runoff_input_mm here and explain why it wouldn’t make sense to use weighted_precip.mm.
ANSWER:

By subtracting discharge and PET from runoff_input, we’re essentially examining the residual water. This could represent the amount of water available to the system after accounting for the demand (PET) and the outflow (discharge).

Q4 (4pnts) Consider the evapotranspiration module and the discussion potential and actual evapotranspiration. Why do you think PET yields estimates of negative residual water in this high elevation, semi-arid system? ANSWER:

Again, you may have to change some of the script below to accommodate your summary dataframe, and not every column will be helpful for a synopsis figure, but the objective in this chunk is to present the column values of your summary dataframe as a bar chart.

Single model run

Look at the parameters to make sure you know what each parameter does. The way the model is set up, is that everything is hidden inside a function. The user (–> you) only formats the input data and passes the parameters to the function. All model output is contained in “modelRun”. If you look at the Environment, you will notice that modelRun is a list. A list is even more flexible in terms of data storage than a dataframe (a dataframe is actually a special type of list…). While lists are super flexible, they can also be more cumbersome to deal with. I included some code that takes the output from the model run and saves all the important bits and pieces in a convenient dataframe called HBVRun. For this part (single model runs), you will really only need the data contained in the HBVRun dataframe.

single model execution

Objective functions

KGE

Before we can start trying to tune our model to look more like the observed discharge record, it would be helpful to have some sort of quantified metric for how well our modeled data fits the measured data.

There are many different ways to do this, but discussion of the pros and cons of those approaches is beyond this quick introduction to modeling.

Here we will demonstrate the Kling-Gupta efficiency both for runoff as well as for swe.

NSE

Nash-Sutcliffe Efficiency (NSE).

Basically, the NSE looks at how much better your model run did that if you had just used the mean discharge for the data record as your “modelled results”. It does this by comparing how far off the observed values where from the mean discharge to how far off the modeled values were from the observed discharge.

Mathematically, NSE is the sum of the squared differences between the modeled and observed discharge divided by the sum of the squared differences between the observed and mean discharge, subtracted by 1.

\[ NSE = 1 - \frac{\sum_{t = 1}^{T}{(Q_m^t - Q_o^t)^2}}{\sum_{t = 1}^{T}{(Q_o^t - \bar{Q_o})^2}} \] Where \(Q_m^t\) is modeled discharge at time t, \(Q_o^t\) is observed discharge at time t, and \(\bar{Q_o}\) is mean observed discharge.

Calibrate HBV manually

Woohoo! We can now run our model and assess how well it is working!

Now, let’s see how well we can get it to work. The code below runs the model, produces a plot, and calculates the NSE based on discharge.

# set up the parameter vector
model_params <- c(
  1.05, # SCF snow correction factor [-] (e.g., 0.9-1.5);
  1.80, # DDF degree day factor [mm/degC/timestep] (e.g., 0.0-5.0 mm/degC/day);
  2, # Tr threshold temperature above which precipitation is rain [degC] (e.g., 1.0-3.0 degC);
  0, # Ts threshold temperature below which precipitation is snow [degC] (e.g., -3.0-1.0 degC);
  -0.336, # Tm threshold temperature above which melt starts [degC] (e.g., -2.0-2.0 degC);
  0.2, # LPrat parameter related to the limit for potential evaporation [-] (e.g., 0.0-1.0);
  121, # FC field capacity, i.e., max soil moisture storage [mm] (e.g., 0-600 mm);
  2.52, # BETA the non linear parameter for runoff production [-] (e.g., 0.0-20.0);
  0.473, # k0 storage coefficient for very fast response [timestep] (e.g., 0.0-2.0 days);
  9.06, # k1 storage coefficient for fast response [timestep] (e.g., 2.0-30.0 days);
  142, # k2 storage coefficient for slow response [timestep] (e.g., 30.0-250.0 days);
  50.1, # lsuz threshold storage state, i.e., the very fast response start if exceeded [mm] (e.g., 1.0-100.0 mm);
  2.38, # cperc constant percolation rate [mm/timestep] (e.g., 0.0-8.0 mm/day);
  10, # bmax maximum base at low flows [timestep] (e.g., 0.0-30.0 days);
  25 # croute free scaling parameter [timestep^2/mm] (e.g., 0.0-50.0 days^2/mm);
)


# set time period
model_in <- indata %>%
  filter(date >= as_date("2017-10-01") & date <= as_date("2022-09-30"))

# set up the model
## THIS IS THE ACTUAL MODEL EXECUTION
modelRun <- TUWmodel(
  prec = model_in$weighted_precip.mm, # precip input
  airt = model_in$Tmean_c, # air temp input
  ep = model_in$PET_mm, # pet input
  area = 1, # one zone for the entire watershed
  param = model_params # input model parameters
)


# get all outputs into a nice df
HBVRun <- tibble(
  Date = model_in$Date, # date
  P_mm = modelRun$prec, # precip
  Tair = modelRun$airt, # air temp
  SWEobs = model_in$weighted_swe.mm, # observed swe
  PET = modelRun$ep,# pet
  Qobs = model_in$Q_mm_day, # observed discharge
  Qsim = modelRun$q[1, ], # simulated discharge
  Qsurf = modelRun$q0[1, ], # surface runoff
  Qsubsurf = modelRun$q1[1, ], # subsurface flow
  Qbase = modelRun$q2[1, ], # groundwater flow
  Rain = modelRun$rain[1, ], # simulated rain
  Snow = modelRun$snow[1, ], # simulated snowfall
  Melt = modelRun$melt[1, ], # simulated melt
  SWEsim = modelRun$swe[1, ], # simulated swe
  Soilmoist = modelRun$moist[1, ], # simulated soil storage
  AET = modelRun$eta[1, ], # simulated evapotranspiration
  StorageUpper = modelRun$suz[1, ], # upper storage value
  StorageLower = modelRun$slz[1, ] # lower storage value
)

#Trim out the warm up period
OutTrim <- HBVRun %>% slice(366:n())

#Calculate NSE
NSE <- 1 - ((sum((OutTrim$Qsim - OutTrim$Qobs) ^ 2)) / 
                 sum((OutTrim$Qobs - mean(OutTrim$Qobs)) ^ 2))

print(NSE)

[1] -0.3496303

Plot observed and model simulations

Generate plots that include and compare the different modeled fluxes from your best NSE. Some of those fluxes can be immediately compared to observed data (e.g., runoff or SWE), while others only exist in simulated form (e.g., storages or outflows of the various runoff components) and need to be assessed with the perceptual model in mind.
Make sure that the axes are properly labeled when you create plots. The script below will get you started.

# Add date to the HBVRun result tibble
HBVRun <- as.data.frame(HBVRun)
HBVRun$Date <- model_in$date
str(HBVRun)

'data.frame':   1826 obs. of  18 variables:
 $ P_mm        : num  26.78 7.65 1.91 0 0 ...
 $ Tair        : num  0.8 -1.6 -0.7 6.2 6.65 0.25 3.2 2.5 -7.8 -1.8 ...
 $ SWEobs      : num  13.3 46.7 53.4 53.4 51.2 ...
 $ PET         : num  1.3 1.15 1.23 1.84 1.76 ...
 $ Qobs        : num  0.282 0.243 0.234 0.23 0.231 ...
 $ Qsim        : num  0.0135 0.0418 0.0661 0.0531 0.0398 ...
 $ Qsurf       : num  0 0 0 0 0 0 0 0 0 0 ...
 $ Qsubsurf    : num  0.1343 0 0 0 0.0421 ...
 $ Qbase       : num  0.0342 0.0434 0.0431 0.0569 0.0732 ...
 $ Rain        : num  10.7 0 0 0 0 ...
 $ Snow        : num  16.07 7.65 1.91 0 0 ...
 $ Melt        : num  2.04 0 0 11.76 12.57 ...
 $ SWEsim      : num  14.825 22.858 24.866 13.101 0.526 ...
 $ Soilmoist   : num  60.1 60.1 60.1 68 75.9 ...
 $ AET         : num  1.3 0 0 1.84 1.76 ...
 $ StorageUpper: num  1.36 0 0 0 0.52 ...
 $ StorageLower: num  4.85 6.16 6.12 8.08 10.39 ...
 $ Date        : Date, format: "2017-10-01" ...

# Round the NSE for display
nse_label <- paste("NSE =", round(NSE, 3))

q_plot <- ggplot(data = HBVRun) +
  geom_line(aes(x = Date, y = Qobs, color = "Qobs")) +  # Qobs
  geom_line(aes(x = Date, y = Qsim, color = "Qsim")) +  # Qsim
  annotate("text",
           x = max(HBVRun$Date) - 100,  # Move label # days from the end
           y = max(HBVRun$Qobs, na.rm = TRUE),
           label = nse_label,
           hjust = 1, vjust = 1.5, size = 5, fontface = "bold") +
  labs(x = NULL, y = "Q (mm/day)", color = NULL, title = "Qobs and Qsim")


# Make it interactive with plotly
ggplotly(q_plot)

NA

Q6. (5 pts) What is the role of the model parameters in representing hydrological processes? What is one effective method to calibrate these parameters without relying on manual adjustment?

Answer:

More exploratory/demonstrative plots

# Swe
swe_plot <- ggplot(data = HBVRun) +
  geom_line(aes(x = Date, y = SWEobs, color = "SWEobs")) + # SWEobs
  geom_line(aes(BLANK)) + # SWEsim
  labs(x = {}, y = "SWE (mm)", color = {})
ggplotly(swe_plot)

Error: object 'BLANK' not found

# Q0, Q1, Q2
q_bucket_plot <- ggplot(data = HBVRun) +
  geom_line(aes(x = Date, y = Qsurf, color = "Qsurf")) + # Qsurf
  geom_line(aes(x = Date, y = Qsubsurf, color = "Qsubsurf")) + #Qsubsurf
  geom_line(aes(x = Date, y = Qbase, color = "Qbase")) + # Qbase
  labs(x = {}, y = "Q (mm)", color = {}, title = "Q0, Q1, and Q2")
ggplotly(q_bucket_plot)

# PET and AET
pet_plot <- ggplot(data = HBVRun) +
  geom_line(aes(x = Date, y = PET, color = "PET")) + # PET
  geom_line(aes(x = Date, y = AET, color = "AET")) + # AET
  labs(x = {}, y = "Flux (mm)", color = {}, title = "PET and AET")
ggplotly(pet_plot)

# Storages
storage_plot <- ggplot(data = HBVRun) +
  geom_line(aes(x = Date, y = Soilmoist, color = "Soil moisture storage")) + # soil moisture storage
  geom_line(aes(x = Date, y = StorageUpper, color = "Upper storage")) + # upper storage
  geom_line(aes(x = Date, y = StorageLower, color = "Lower storage")) + # lower storage
  labs(x = {}, y = "Storage (mm)", color = {}, title = "Storages")
ggplotly(storage_plot)